#include <stdio.h>
#include <random>
#include <cblas.h>

#define AMAX BLASFUNC(idamax)
#define AMIN BLASFUNC(idamin)
#define DOT BLASFUNC(ddot)
#define NRM2 BLASFUNC(dnrm2)
#define ASUM BLASFUNC(dasum)
#define SCAL BLASFUNC(dscal)
#define AXPY BLASFUNC(daxpy)

#define FLOAT double
#define SCALAR double
#define SCALAR_SUFFIX d
#define EIGEN_NO_DEBUG
#define EIGEN_DONT_PARALLELIZE

#include "blas_level1_impl.h"
#include "blas_level1_real_impl.h"

struct timespec start1, stop1;
struct timespec start2, stop2;

void print_time(const char *title, double total1, double total2, int loops)
{
    fprintf(stderr, title);
    fprintf(stderr,"Openblas average time: %10.9f s\n", total1 / loops);
    fprintf(stderr, "Eigen 3 average time: %10.9f s\n", total2 / loops);
    fprintf(stderr,"diff = %f%%\n", 100 * (total2 - total1) / total1);
}

void amax(const int n, const int &loops)
{
    FLOAT *x = new FLOAT[n];
    double total1 = 0.0;
    double total2 = 0.0;
    int inc_x = 1;

    for (size_t i = 0; i < loops; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            x[j] = ((FLOAT)rand() / (FLOAT)RAND_MAX) - 0.5;
        }
        clock_gettime(CLOCK_REALTIME, &start1);
        int result1 = AMAX((int *)&n, x, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop1);

        clock_gettime(CLOCK_REALTIME, &start2);
        int result2 = eigen_amax((int *)&n, x, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop2);

        double time1 = (double)(stop1.tv_sec - start1.tv_sec) + (double)((stop1.tv_nsec - start1.tv_nsec)) * 1.e-9;
        total1 += time1;

        double time2 = (double)(stop2.tv_sec - start2.tv_sec) + (double)((stop2.tv_nsec - start2.tv_nsec)) * 1.e-9;
        total2 += time2;

        if(result1!=result2)
        {
            fprintf(stderr, "AMAX FAIL %d %d\n", result1, result2);
        }
    }
    print_time("---AMAX---\n", total1, total2, loops);

    delete[] x;
}

void amin(const int n, const int &loops)
{
    FLOAT *x = new FLOAT[n];
    FLOAT total1 = 0.0;
    FLOAT total2 = 0.0;
    int inc_x = 1;

    for (size_t i = 0; i < loops; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            x[j] = ((FLOAT)rand() / (FLOAT)RAND_MAX) - 0.5;
        }
        clock_gettime(CLOCK_REALTIME, &start1);
        int result1 = AMIN((int *)&n, x, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop1);

        clock_gettime(CLOCK_REALTIME, &start2);
        int result2 = eigen_amin((int *)&n, x, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop2);

        double time1 = (double)(stop1.tv_sec - start1.tv_sec) + (double)((stop1.tv_nsec - start1.tv_nsec)) * 1.e-9;
        total1 += time1;

        double time2 = (double)(stop2.tv_sec - start2.tv_sec) + (double)((stop2.tv_nsec - start2.tv_nsec)) * 1.e-9;
        total2 += time2;

        if(result1!=result2)
        {
            fprintf(stderr, "AMIN FAIL  %d %d\n", result1, result2);
        }
    }

    print_time("---AMIN---\n", total1, total2, loops);

    delete[] x;
}

void dot(const int n, const int &loops)
{
    FLOAT *x = new FLOAT[n];
    FLOAT *y = new FLOAT[n];
    FLOAT total1 = 0.0;
    FLOAT total2 = 0.0;
    int inc_x = 1;

    for (size_t i = 0; i < loops; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            x[j] = ((FLOAT)rand() / (FLOAT)RAND_MAX) - 0.5;
            y[j] = ((FLOAT)rand() / (FLOAT)RAND_MAX) - 0.5;
        }
        clock_gettime(CLOCK_REALTIME, &start1);
        FLOAT result1 = DOT((int *)&n, x, &inc_x, y, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop1);

        clock_gettime(CLOCK_REALTIME, &start2);
        FLOAT result2 = eigen_dot((int *)&n, x, &inc_x, y, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop2);

        double time1 = (double)(stop1.tv_sec - start1.tv_sec) + (double)((stop1.tv_nsec - start1.tv_nsec)) * 1.e-9;
        total1 += time1;

        double time2 = (double)(stop2.tv_sec - start2.tv_sec) + (double)((stop2.tv_nsec - start2.tv_nsec)) * 1.e-9;
        total2 += time2;

        if(abs(result1-result2)>1e-8)
        {
            fprintf(stderr, "DOT FAIL %f %f\n", result1, result2);
        }
    }

    print_time("---DOT---\n", total1, total2, loops);

    delete[] x;
    delete[] y;
}

void nrm2(const int n, const int &loops)
{
    FLOAT *x = new FLOAT[n];
    FLOAT total1 = 0.0;
    FLOAT total2 = 0.0;
    int inc_x = 1;

    for (size_t i = 0; i < loops; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            x[j] = ((FLOAT)rand() / (FLOAT)RAND_MAX) - 0.5;
        }
        clock_gettime(CLOCK_REALTIME, &start1);
        FLOAT result1 = NRM2((int *)&n, x, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop1);

        clock_gettime(CLOCK_REALTIME, &start2);
        FLOAT result2 = eigen_nrm2((int *)&n, x, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop2);

        double time1 = (double)(stop1.tv_sec - start1.tv_sec) + (double)((stop1.tv_nsec - start1.tv_nsec)) * 1.e-9;
        total1 += time1;

        double time2 = (double)(stop2.tv_sec - start2.tv_sec) + (double)((stop2.tv_nsec - start2.tv_nsec)) * 1.e-9;
        total2 += time2;

        if(abs(result1-result2)>1e-8)
        {
            fprintf(stderr, "NRM2 FAIL  %f %f\n", result1, result2);
        }
    }

    print_time("---NRM2---\n", total1, total2, loops);

    delete[] x;
}

void asum(const int n, const int &loops)
{
    FLOAT *x = new FLOAT[n];
    FLOAT total1 = 0.0;
    FLOAT total2 = 0.0;
    int inc_x = 1;

    for (size_t i = 0; i < loops; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            x[j] = ((FLOAT)rand() / (FLOAT)RAND_MAX) - 0.5;
        }
        clock_gettime(CLOCK_REALTIME, &start1);
        FLOAT result1 = ASUM((int *)&n, x, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop1);

        clock_gettime(CLOCK_REALTIME, &start2);
        FLOAT result2 = eigen_asum((int *)&n, x, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop2);

        double time1 = (double)(stop1.tv_sec - start1.tv_sec) + (double)((stop1.tv_nsec - start1.tv_nsec)) * 1.e-9;
        total1 += time1;

        double time2 = (double)(stop2.tv_sec - start2.tv_sec) + (double)((stop2.tv_nsec - start2.tv_nsec)) * 1.e-9;
        total2 += time2;

        if(abs(result1-result2)>1e-8)
        {
            fprintf(stderr, "ASUM FAIL  %f %f\n", result1, result2);
        }
    }

    print_time("---ASUM---\n", total1, total2, loops);

    delete[] x;
}


void scal(const int n, const int &loops)
{
    FLOAT *x = new FLOAT[n];
    FLOAT total1 = 0.0;
    FLOAT total2 = 0.0;
    int inc_x = 1;
    FLOAT alpha[2] = {2.0, 2.0};

    for (size_t i = 0; i < loops; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            x[j] = ((FLOAT)rand() / (FLOAT)RAND_MAX) - 0.5;
        }
        clock_gettime(CLOCK_REALTIME, &start1);
        SCAL((int *)&n, alpha, x, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop1);

        clock_gettime(CLOCK_REALTIME, &start2);
        eigen_scal((int *)&n, alpha, x, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop2);

        double time1 = (double)(stop1.tv_sec - start1.tv_sec) + (double)((stop1.tv_nsec - start1.tv_nsec)) * 1.e-9;
        total1 += time1;

        double time2 = (double)(stop2.tv_sec - start2.tv_sec) + (double)((stop2.tv_nsec - start2.tv_nsec)) * 1.e-9;
        total2 += time2;
    }

    print_time("---SCAL---\n", total1, total2, loops);

    delete[] x;
}

void axpy(const int n, const int &loops)
{
    FLOAT *x = new FLOAT[n];
    FLOAT *y = new FLOAT[n];
    FLOAT total1 = 0.0;
    FLOAT total2 = 0.0;
    int inc_x = 1;
    FLOAT alpha[2] = {2.0, 2.0};

    for (size_t i = 0; i < loops; i++)
    {
        for (size_t j = 0; j < n; j++)
        {
            x[j] = ((FLOAT)rand() / (FLOAT)RAND_MAX) - 0.5;
            y[j] = ((FLOAT)rand() / (FLOAT)RAND_MAX) - 0.5;
        }
        clock_gettime(CLOCK_REALTIME, &start1);
        AXPY((int *)&n, alpha, x, &inc_x, y, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop1);

        clock_gettime(CLOCK_REALTIME, &start2);
        eigen_axpy((int *)&n, alpha, x, &inc_x, y, &inc_x);
        clock_gettime(CLOCK_REALTIME, &stop2);

        double time1 = (double)(stop1.tv_sec - start1.tv_sec) + (double)((stop1.tv_nsec - start1.tv_nsec)) * 1.e-9;
        total1 += time1;

        double time2 = (double)(stop2.tv_sec - start2.tv_sec) + (double)((stop2.tv_nsec - start2.tv_nsec)) * 1.e-9;
        total2 += time2;
    }

    print_time("---AXPY---\n", total1, total2, loops);

    delete[] x;
}


int main(int argc, char **argv)
{
    int n = std::stoi(argv[1]);
    int loops = std::stoi(argv[2]);

    #ifdef linux
    srand(getpid());
    #endif

    amax(n, loops);
    amin(n, loops);
    dot(n, loops);
    nrm2(n, loops);
    asum(n, loops);
    scal(n, loops);
    axpy(n, loops);

    return 0;
}